1. Ejercicio 1¶
library("IRdisplay")
# Enseño el grafo dado por el experto
display_png(file="grafo.png")
- A partir del grafo dado por el experto, determina la expresión que refleja la factorización de la distribución global, y escríbela como una secuencia de caracteres del tipo generado por la función modelstring del paquete de R bnlearn.
NOTA: Como ya se comentado anteriormente, es importante evitar tildes a la hora de nombrar los nodos. Tampoco son recomendables espacios en blanco ni ningún tipo de carácter especial.
A continuación, introduce el DAG en R utilizando la definicíon del modelo que acabas de crear
Realiza una lista de los padres e hijos de cada uno de los nodos
Realiza una lista de todas las conexiones fundamentales presentes en el grafo, y realiza una clasificación de cada una de ellas determinando si es una estructura en serie, divergente o convergente.
¿Se observa alguna v-estructura en el grafo?
Intenta introducir un arco que parta del nodo neblina y se dirija hacia escarcha, y otro que vaya de granizo a lluvia. Comenta qué sucede en cada caso, y si el resultado daría lugar a una red bayesiana válida.
Determina la manta de Markov del nodo rocio
Introduce un arco que parta del nodo Lluvia y se dirija al nodo Niebla, ¿cambia la manta de Markov del nodo rocio? En caso afirmativo, ¿cómo lo hace?
Para continuar, elimina el arco definido anteriormente entre los nodos Lluvia y Niebla, recuperando el DAG original.
Desarrollo ejercicio 1¶
A partir del grafo dado por el experto, determina la expresión que refleja la factorización de la distribución global, y escríbela como una secuencia de caracteres del tipo generado por la función modelstring del paquete de R bnlearn.¶
[Viento][Rocio][Escarcha][Niebla][Nieblina][Lluvia][Granizo][Tormenta][Nieve][Nieve_Suelo]
A continuación, introduce el DAG en R utilizando la definicíon del modelo que acabas de crear¶
library(bnlearn)
## Defino el grafo vacío:
dag <- empty.graph(nodes=c("viento", "rocio", "escarcha", "niebla", "neblina",
"lluvia", "granizo", "tormenta", "nieve", "nieveSuelo"))
# Defino los arcos
arc.set <- matrix(c("viento", "lluvia",
"viento", "niebla",
"viento", "rocio",
"rocio", "niebla",
"rocio", "escarcha",
"niebla", "neblina",
"lluvia", "tormenta",
"tormenta", "granizo",
"tormenta", "nieve",
"nieve", "nieveSuelo"),
byrow = TRUE, ncol = 2,
dimnames = list(NULL, c("from", "to")))
arcs(dag) <- arc.set
print(dag)
plot(dag)
# Verifico la factorización del punto anterior
modelstring(dag)
Realiza una lista de los padres e hijos de cada uno de los nodos¶
for (node in nodes(dag)){
print(paste("NODO: ", node))
print("----")
print("PADRES:")
print(parents(dag, node = node))
print("----")
print("HIJOS:")
print(children(dag, node = node))
print("====")
}
Realiza una lista de todas las conexiones fundamentales presentes en el grafo, y realiza una clasificación de cada una de ellas determinando si es una estructura en serie, divergente o convergente.¶
- Estructuras convergentes:
Se encuentra la estructura convergente
viento --> niebla <-- rocio
Por otro lado, esa no representa una v-estrucura, porque hay un arco que conecta viento con rocio.
- Estructuras divergentes:
El grafo presenta 3 estructuras divergentes.
- De viento salen 3 arcos:
- viento --> rocio
- viento --> niebla
- viento --> lluvia
- De rocio salen 2 arcos:
- rocio --> niebla
- rocio --> escarcha
- De tormenta salen 2 arcos:
- tormenta --> nieve
- tormenta --> granizo
Estructuras en serie:
- viento --> niebla --> nieblina
- viento --> rocio --> niebla --> nieblina
- viento --> rocio --> escarcha
- viento --> lluvia --> tormenta --> nieve --> nieveSuelo
- viento --> lluvia --> tormenta --> granizo
¿Se observa alguna v-estructura en el grafo?¶
No hay v-estructuras en el grafo
# Verifico que no haya v-estructuras
vstructs(dag)
Intenta introducir un arco que parta del nodo neblina y se dirija hacia escarcha, y otro que vaya de granizo a lluvia. Comenta qué sucede en cada caso, y si el resultado daría lugar a una red bayesiana válida.¶
- Un arco que conecta el nodo "neblina" con el nodo "escarcha" introduce una nueva v-estructura:
neblina --> escarcha <-- rocío
dag <- set.arc(dag, from = "neblina", to = "escarcha")
vstructs(dag)
plot(dag)
- Un arco que conecte el nodo "Granizo" con el nodo "Lluvia" crearía un ciclo, y por eso está prohibo. Vemos qué dice el código si intentamos introducirlo.
dag <- set.arc(dag, from = "granizo", to = "lluvia")
Determina la manta de Markov del nodo rocio¶
mb(dag, "rocio")
Introduce un arco que parta del nodo Lluvia y se dirija al nodo Niebla, ¿cambia la manta de Markov del nodo rocio? En caso afirmativo, ¿cómo lo hace?¶
La introducción de este arco genera una nueva v-estrucura: rocio --> niebla <-- lluvia.
Esto hace que el nodo lluvia entre en la manta de Markov del nodo rocio.
dag <- set.arc(dag, from = "lluvia", to = "niebla")
vstructs(dag)
mb(dag, "rocio")
plot(dag)
Para continuar, elimina el arco definido anteriormente entre los nodos Lluvia y Niebla, recuperando el DAG original.¶
# Vuelvo a definir los arcos originales
arc.set <- matrix(c("viento", "lluvia",
"viento", "niebla",
"viento", "rocio",
"rocio", "niebla",
"rocio", "escarcha",
"niebla", "neblina",
"lluvia", "tormenta",
"tormenta", "granizo",
"tormenta", "nieve",
"nieve", "nieveSuelo"),
byrow = TRUE, ncol = 2,
dimnames = list(NULL, c("from", "to")))
arcs(dag) <- arc.set
print(dag)
plot(dag)
2. Ejercicio 2¶
Considerando la información proporcionada por el conjunto de datos meteor, construye la red bayesiana utilizando el método bayesiano de estimación paramétrica.
¿Cuál sería el número potencial de parámetros (sin usar la red bayesiana) del modelo para calcular la probabilidad global si no utilizasemos el DAG?
¿Cuántos parámetros tiene la distribución global dada por la red bayesiana?
Determina el número de parámetros asociado a cada una de las distribuciones locales
Obtén las tablas de probabilidad condicional asociadas los nodos granizo y niebla. Ahora representa la información de cada tabla en sendos gráficos.
Desarrollo ejercicio 2¶
Considerando la información proporcionada por el conjunto de datos meteor, construye la red bayesiana utilizando el método bayesiano de estimación paramétrica.¶
# Construyo la red bayesiana a partir de la tabla
meteoro <- read.table("meteoro.txt", header = TRUE)
bn = bn.fit(dag, data = meteoro, method = "bayes")
¿Cuál sería el número potencial de parámetros (sin usar la red bayesiana) del modelo para calcular la probabilidad global si no utilizasemos el DAG?¶
número de parametros = p = $2^n$ - 1, con n = numero de nodos = 10
p = $2^{10} - 1$ = 1023
Considerando el grafo, que hemos contruido, los parametros son dados por la factorización:
[viento][rocio|viento][lluvia|viento][escarcha|rocio][niebla|viento:rocio] [tormenta|lluvia][neblina|niebla][granizo|tormenta][nieve|tormenta] [nieveSuelo|nieve]
Miramos cuantos parametros tiene cada uno de los factores:
[viento]: 1 parametro (probabilidad de 's' y probabilidad de 'n' = 1 - p('s'))
[rocio|viento]: 2 parametros
[lluvia|viento]: 2 parametros
[escarcha|rocio]: 2 parametros
[niebla|viento:rocio] : 4 parametros
[tormenta|lluvia]: 2 parametros
[neblina|niebla]: 2 parametros
[granizo|tormenta]: 2 parametros
[nieve|tormenta]: 2 parametros
[nieveSuelo|nieve]: 2 parametros
Sumando: 21
# Verifico con el código
nparams(bn)
Obtén las tablas de probabilidad condicional asociadas a los nodos granizo y niebla. Ahora representa la información de cada tabla en sendos gráficos.¶
# Nodo granizo
print(bn$granizo$prob)
plot(bn$granizo$prob)
# Nodo niebla
print(bn$niebla$prob)
plot(bn$niebla$prob)
3. Ejercicio 3¶
Una vez construída la red hemos establecido la base de conocimiento del sistema inteligente. A continuación se puede calcular la probabilidad de cualquier variable o conjunto de variables condicionadas a cualquier evidencia que se tenga disponible para un problema dado, es decir, realizar la inferencia.
- Observa la estructura del DAG y repasa la teoría y el concepto de d-separación. Responde razonadamente si las siguientes afirmaciones son verdaderas o falsas, utilizando únicamente la estructura del DAG y el concepto de d-separación:
- La nieve y el granizo son fenómenos independientes a priori
- La nieve y el granizo son fenómenos independientes dado que haya habido tormenta
- La nieve en el suelo y la neblina son fenómenos independientes
- La nieve en el suelo y la neblina son fenómenos independientes dado que haya habido tormenta
NOTA: En este apartado se deberá aplicar la inferencia exacta. Conocido que en un día dado se han producido tormentas, calcula cómo afecta este hecho a la probabilidad de que se produzcan los siguientes fenómenos meteorológicos:
- Que llueva, P(lluvia=s|tormenta=s)
- Que haya rachas de viento superiores a 50 Km/h, P(viento=s|tormenta=s)
- Que llueva y que además las rachas de viento superen los 50 Km/h, P(lluvia=s,viento=s|tormenta=s)
- A partir de la información revelada por la red bayesiana, sabiendo que un día se producen tormentas:
- ¿Hay mayor probabilidad de que llueva cuando se producen tormentas que cualquier otro día?
- ¿Aumenta o disminuye la probabilidad de tener rachas de viento mayores de 50 Km/h cuando se produce tormenta? ¿Cuánto?
Dado el ejercicio anterior, repetir ahora el ejercicio mediante inferencia aproximada, calculando para cada una de las estimaciones 100 realizaciones. Para cada una de las respuestas anteriores, representa un diagrama de cajas que represente la dispersión en la estimación de la probabilidad, marcando además el valor obtenido de manera exacta en el apartado anterior.
Desarrollo ejercicio 3¶
# Vuelvo a pintar el grafo
plot(dag)
Observa la estructura del DAG y repasa la teoría y el concepto de d-separación. Responde razonadamente si las siguientes afirmaciones son verdaderas o falsas, utilizando únicamente la estructura del DAG y el concepto de d-separación:¶
La nieve y el granizo son fenómenos independientes a priori
No, ya que tienen el mismo padre y no son d-separados. Podrían ser independientes, pero no lo podemos afirmar con lo que conocemos.
# Verifico
dsep(dag, x = "nieve", y = "granizo")
La nieve y el granizo son fenómenos independientes dado que haya habido tormenta
Sí, en este caso el conocer el estado del nodo tormenta hace que el camino entre nieve y tormenta no esté activo. Por lo tanto, los dos nodos son d-separados y entonces independientes.
# Verifico
dsep(dag, x = "nieve", y = "granizo", z = "tormenta")
La nieve en el suelo y la neblina son fenómenos independientes
No: como existen dos caminos que las conectan, podrían ser independientes, pero no lo podemos afirmar con lo que conocemos.
# Verifico
dsep(dag, x = "nieveSuelo", y = "neblina")
La nieve en el suelo y la neblina son fenómenos independientes dado que haya habido tormenta
Sí, como en caso de nieve y granizo, si conocemos el estado del nodo tormenta se interrumpe el camino entre niebla y nieve en el suelo, y las dos variables se d-separan.
# Verifico
dsep(dag, x = "nieveSuelo", y = "neblina", z = "tormenta")
NOTA: En este apartado se deberá aplicar la inferencia exacta. Conocido que en un día dado se han producido tormentas, calcula cómo afecta este hecho a la probabilidad de que se produzcan los siguientes fenómenos meteorológicos:¶
library(gRain)
# Compilo la red
red <- compile(as.grain(bn))
# Fijo el estado del nodo tormenta a 's'
red.tormenta.s = setEvidence(red, nodes = "tormenta", states = "s")
- Que llueva, P(lluvia=s|tormenta=s)
# Consulto la red: P(lluvia=s|tormenta=s) = ?
query1 <- querygrain(red.tormenta.s, nodes = "lluvia", type = "marginal")$lluvia["s"]
print(paste("La probabilidad de lluvia dado que hay tormenta es", query1))
- Que haya rachas de viento superiores a 50 Km/h, P(viento=s|tormenta=s)
# Consulto la red: P(viento=s|tormenta=s) = ?
query2 <- querygrain(red.tormenta.s, nodes = "viento", type = "marginal")$viento["s"]
print(paste("La probabilidad de viento dado que hay tormenta es", query2))
- Que llueva y que además las rachas de viento superen los 50 Km/h, P(lluvia=s,viento=s|tormenta=s)
# Consulto la red: P(lluvia=s,viento=n|tormenta=s) = ?
query3 <- querygrain(red.tormenta.s, nodes = c("lluvia","viento"), type = "joint")["s","s"]
print(paste("La probabilidad de viento dado que hay tormenta es", query3))
- A partir de la información revelada por la red bayesiana, sabiendo que un día se producen tormentas:
- ¿Hay mayor probabilidad de que llueva cuando se producen tormentas que cualquier otro día?
# Miro la probabilidad de que llueva, sin poner condiciones
query4 <- querygrain(red, nodes = "lluvia", type = "marginal")$lluvia["s"]
print(paste("La probabilidad de lluvia dado que hay tormenta es", query1,
"mientras que sin poner condiciones, la probabilidad de que llueva es", query4,
". Eso quiere decir que es más probable que llueva cuando hay tormenta."))
- Y ahora:
- ¿Aumenta o disminuye la probabilidad de tener rachas de viento mayores de 50 Km/h cuando se produce tormenta? ¿Cuánto?
# Miro la probabilidad de que haya viento, sin poner condiciones
query5 <- querygrain(red, nodes = "viento", type = "marginal")$viento["s"]
print(paste("La probabilidad de que haya viento dado que hay tormenta es", query2,
"mientras que sin poner condiciones, la probabilidad de que haya viento es", query5,
". Eso quiere decir que es más probable que haya viento cuando hay tormenta.",
"Cuantitativamente, el aumento es del", 100*(query2 / query5 - 1), "%"))
- Dado el ejercicio anterior, repetir ahora el ejercicio mediante inferencia aproximada, calculando para cada una de las estimaciones 100 realizaciones. Para cada una de las respuestas anteriores, representa un diagrama de cajas que represente la dispersión en la estimación de la probabilidad, marcando además el valor obtenido de manera exacta en el apartado anterior.
# Fijo el numero de simulaciones a 100
N = 100
# Simulo 100 veces lo que pasa a la red sin poner condiciones
bn.sim <- simulate(red, n = N)
# Verifico que tormenta tenga como valores tantos 's' como 'n'
print(paste("Número de 's'", nrow(bn.sim[bn.sim$tormenta == 's',])))
print(paste("Número de 'n'", nrow(bn.sim[bn.sim$tormenta == 'n',])))
# Miro la probabilidad de que haya tormenta, sin poner condiciones
# para verificar que los números que saco tienen sentido
query.tormenta <- querygrain(red, nodes = "tormenta", type = "marginal")$tormenta
query.tormenta
# Simulo 100 veces lo que pasa a la red si siempre hay tormenta
bn.sim.tormenta.s <- simulate(red.tormenta.s, n = 100)
# Verifico que tormenta siempre tenga como valor 's'
bn.sim.tormenta.s$tormenta
Ahora que tengo las simulaciones, puedo hacer consultas:
- Que llueva, P(lluvia=s|tormenta=s)
prob.lluvia <- nrow(bn.sim.tormenta.s[bn.sim.tormenta.s$lluvia == 's',])/N
print(paste("Probabilidad de lluvia, dada tormenta:",
prob.lluvia))
my.table = table(bn.sim.tormenta.s$lluvia)
plot(my.table, type = "p", ylim = c(0,100))
arrows(1, my.table["n"] + sqrt(my.table["n"]),
1, my.table["n"] - sqrt(my.table["n"]),
length=0.05, angle=90, code=3)
points(1, 100*(1 - query1), col = "red")
arrows(2, my.table["s"] + sqrt(my.table["s"]),
2, my.table["s"] - sqrt(my.table["s"]),
length=0.05, angle=90, code=3)
points(2, 100*query1, col = "red")
- Que haya rachas de viento superiores a 50 Km/h, P(viento=s|tormenta=s)
prob.viento <- nrow(bn.sim.tormenta.s[bn.sim.tormenta.s$viento == 's',])/N
print(paste("Probabilidad de viento, dada tormenta:",
prob.viento))
my.table = table(bn.sim.tormenta.s$viento)
plot(my.table, type = "p", ylim = c(0,100))
arrows(1, my.table["n"] + sqrt(my.table["n"]),
1, my.table["n"] - sqrt(my.table["n"]),
length=0.05, angle=90, code=3)
points(1, 100*(1 - query2), col = "red")
arrows(2, my.table["s"] + sqrt(my.table["s"]),
2, my.table["s"] - sqrt(my.table["s"]),
length=0.05, angle=90, code=3)
points(2, 100*query2, col = "red")
Que llueva y que además las rachas de viento superen los 50 Km/h, P(lluvia=s,viento=s|tormenta=s)
prob.viento.lluvia <- nrow(bn.sim.tormenta.s[bn.sim.tormenta.s$viento == 's' & bn.sim.tormenta.s$lluvia == 's',])/N
print(paste("Probabilidad de viento y lluvia, dada tormenta:",
prob.viento.lluvia))
my.table = table(bn.sim.tormenta.s$viento == 's' & bn.sim.tormenta.s$lluvia == 's')
names(my.table) = c("n", "s")
plot(my.table, type = "p", ylim = c(0,100))
arrows(1, my.table["n"] + sqrt(my.table["n"]),
1, my.table["n"] - sqrt(my.table["n"]),
length=0.05, angle=90, code=3)
points(1, 100*(1 - query3), col = "red")
arrows(2, my.table["s"] + sqrt(my.table["s"]),
2, my.table["s"] - sqrt(my.table["s"]),
length=0.05, angle=90, code=3)
points(2, 100*query3, col = "red")
- A partir de la información revelada por la red bayesiana, sabiendo que un día se producen tormentas:
- ¿Hay mayor probabilidad de que llueva cuando se producen tormentas que cualquier otro día?
print(paste("La probabilidad de lluvia dado que hay tormenta es",
prob.lluvia))
prob.lluvia.no_cond <- nrow(bn.sim[bn.sim$lluvia == 's',])/N
print(paste("Mientras que sin poner condiciones, la probabilidad de que llueva es",
prob.lluvia.no_cond))
print("Eso quiere decir que es más probable que llueva cuando hay tormenta.")
my.table = table(bn.sim$lluvia)
plot(my.table, type = "p", ylim = c(0,100))
arrows(1, my.table["n"] + sqrt(my.table["n"]),
1, my.table["n"] - sqrt(my.table["n"]),
length=0.05, angle=90, code=3)
points(1, 100*(1 - query4), col = "red")
arrows(2, my.table["s"] + sqrt(my.table["s"]),
2, my.table["s"] - sqrt(my.table["s"]),
length=0.05, angle=90, code=3)
points(2, 100*query4, col = "red")
- Y ahora:
- ¿Aumenta o disminuye la probabilidad de tener rachas de viento mayores de 50 Km/h cuando se produce tormenta? ¿Cuánto?
print(paste("La probabilidad de que haya viento dado que hay tormenta es",
prob.viento))
prob.viento.no_cond <- nrow(bn.sim[bn.sim$viento == 's',])/N
print(paste("Mientras que sin poner condiciones, la probabilidad de que haya viento es",
prob.viento.no_cond))
print("Eso quiere decir que es menos probable que haya viento cuando hay tormenta.")
print(paste("Cuantitativamente, disminuye del",
100 * (prob.viento / prob.viento.no_cond - 1), "%"))
my.table = table(bn.sim$viento)
plot(my.table, type = "p", ylim = c(0,100))
arrows(1, my.table["n"] + sqrt(my.table["n"]),
1, my.table["n"] - sqrt(my.table["n"]),
length=0.05, angle=90, code=3)
points(1, 100*(1 - query5), col = "red")
arrows(2, my.table["s"] + sqrt(my.table["s"]),
2, my.table["s"] - sqrt(my.table["s"]),
length=0.05, angle=90, code=3)
points(2, 100*query5, col = "red")
4. Ejercicio 4¶
Como hemos visto, es posible realizar un aprendizaje automático de la estructura del DAG a partir de los datos, usando algoritmos específicos para ello. Durante las clases hemos visto el ejemplo del algoritmo hill-climbing, aunque como vimos, hay otras posibilidades. Tambien hemos visto que podemos combinar nuestra experiencia y el aprendizaje automático definiendo previamente relaciones entre variables que queremos introducir o descartar en el DAG resultante. Además, se ha explicado que existen scores que sirven como criterio para evaluar la fuerza de la dependendencia entre nodos de la red y comparar la bondad de ajuste del modelo.
Evalúa la significación de los arcos dibujados por el experto en el actual DAG utilizando el estadístico χ2 ¿Hay algún arco no significativo? ¿Cuáles son los tres pares de nodos que presentan un arco de unión más fuerte?
Además del algoritmo hill-climbing, existe otro algoritmo popular de tipo “voraz” denominado Tabu search. En bnlearn se encuentra implementado a traves de la función tabu, y los argumentos de entrada son similares a los vistos para hill-climbing.
- Compara el score global (BIC) obtenido por el DAG inicial, con los obtenidos por los DAG aprendidos de forma automática con los algoritmos tabu y hill-climbing. ¿Cuál obtiene mejor score?
- Ahora, para comparar el DAG original con los dos nuevos DAG tabu y hill-climbing, dibuja los tres utilizando la función graphviz.plot. A la luz de los grafos obtenidos en cada caso, ¿cuál te parece que recoge mejor las relaciones causa-efecto entre variables?. Hay que tener en cuenta que los arcos de un grafo no expresan causalidad, sino simplemente dependencia entre variables en términos de probabilidad.
Ahora vuelve a aprender de forma automática el DAG usando tabu y hill-climbing, pero imponiendo las siguientes restricciones: 1. Los arcos viento --> lluvia, tormenta --> granizo y nieve --> nieveSuelo deben quedar reflejados en el DAG. 2. Ningún arco debe unir directamente la neblina con el granizo ni la niebla con la tormenta.
- Vuelve a dibujar los DAG resultantes, y a partir del BIC obtenido por cada modelo, determina cuál es el mejor de todos poniéndolos en una tabla.
- Computa la fuerza de la relación entre nodos de los dos últimos modelos que combinan nuestra experiencia con el aprendizaje automático.
- Comenta brevemente los resultados obtenidos tras combinar nuestro conocimiento con el aprendizaje automático.
- Vuelve a dibujar los DAG resultantes, y a partir del BIC obtenido por cada modelo, determina cuál es el mejor de todos poniéndolos en una tabla.
Desarrollo ejercicio 4¶
- Evalúa la significación de los arcos dibujados por el experto en el actual DAG utilizando el estadístico χ2 ¿Hay algún arco no significativo? ¿Cuáles son los tres pares de nodos que presentan un arco de unión más fuerte?
# Uso arc.strength para valorar la significación de los arcos
arc.strength(dag, data = meteoro, criterion = "x2")
Se observa que los p-valores de todos los arcos son muy pequeños ($< 10^{-5}$) y por eso podemos considerar que todos son significativos. Dicho de otra manera, los arcos sugeridos por el experto tienen en cuenta relaciones probabilisticas que de verdad existen en el conjunto de datos.
En particular, los arcos más significativo (los que presentan p-valor más pequeño) son:
- tormenta --> granizo
- nieve --> nieveSuelo
- viento --> lluvia
- Además del algoritmo hill-climbing, existe otro algoritmo popular de tipo “voraz” denominado Tabu search. En bnlearn se encuentra implementado a traves de la función tabu, y los argumentos de entrada son similares a los vistos para hill-climbing.
# Construyo el DAG usando el algoritmo Tabu search
dag.tabu <- tabu(x = meteoro)
dag.tabu
plot(dag.tabu)
# Construyo el DAG usando el algoritmo hill-climbing
dag.hc <- hc(x = meteoro)
dag.hc
plot(dag.hc)
- Compara el score global (BIC) obtenido por el DAG inicial, con los obtenidos por los DAG aprendidos de forma automática con los algoritmos tabu y hill-climbing. ¿Cuál obtiene mejor score?
# BIC score para el DAG original
dag.score <- bnlearn::score(dag, data = meteoro)
dag.score
# BIC score para el DAG generado con el algoritmo Tabu search
dag.tabu.score <- bnlearn::score(dag.tabu, data = meteoro)
dag.tabu.score
# BIC score para el DAG generado con el algoritmo hill-climbing
dag.hc.score <- bnlearn::score(dag.hc, data = meteoro)
dag.hc.score
El mejor BIC (más cerca de 0, o más pequeño en valor absoluto) se obtiene con el algoritmo Tabu search.
- Ahora, para comparar el DAG original con los dos nuevos DAG tabu y hill-climbing, dibuja los tres utilizando la función graphviz.plot. A la luz de los grafos obtenidos en cada caso, ¿cuál te parece que recoge mejor las relaciones causa-efecto entre variables?. Hay que tener en cuenta que los arcos de un grafo no expresan causalidad, sino simplemente dependencia entre variables en términos de probabilidad.
# graphviz.plot me dá el siguiente problema en local:
graphviz.plot(dag)
# Por esta razón, he dibujado los grafos en el datasciencehub (donde no logro instalar gRain)
# y los pego aquí abajo
print("DAG experto")
display_png(file="DAG.png")
print("DAG tabu")
display_png(file="DAG_Tabu.png")
# Uso arc.strength para valorar la significación de los arcos
arc.strength(dag.tabu, data = meteoro, criterion = "x2")
print("DAG hill-climbing")
display_png(file="DAG_hc.png")
# Uso arc.strength para valorar la significación de los arcos
arc.strength(dag.hc, data = meteoro, criterion = "x2")
Sin duda el grafo más legible e interpretable es el que ha proporcionado el experto, donde se notan relaciones que podemos considerar 'lógicas' y facilmente explicables, como:
- nieve --> nieveSuelo;
- rocio --> escarcha;
- niebla --> neblina.
Por otro lado, este conocimiento, aunque importante y fundamental, puede resultar limitado a la sola experiencia y capacidad de observación de quién estudia el problema, mientras que relaciones importantes entre las caracteristicas de los datos pueden estar escondidas y aparecer solo gracias al utilizo de algoritmos de aprendizaje.
Se observa por ejemplo, que en los dos grafos construidos con algoritmos de aprendizaje, aparece un arco (en un caso muy fuerte) entre nieve y granizo, que no aparece en el grafo experto.
Hay que considerar que la otra cara de este tipo de aprendizaje basado solo en los datos puede considerar relaciones probabilisticas debiles y probablemente basadas más en fluctuaciones estadisticas que en reales dependencia entre dos eventos. Por ejemplo, me resulta extraño que haya relación entre granizo rocio. Se pueden generar así arcos que en vez de describir relaciones entre las caracteristicas del problema, los hacen sobre caracteristicas del dataset (como un sobreajuste), haciendo que el número de parametros del modelo crezca sin aportar realmente nada a su capacidad predictiva. Para obviar a este inconveniente, se suelen introducir términos de regularización en la computación del score del modelo, para penalizar modelos demasiado complejos.
- Ahora vuelve a aprender de forma automática el DAG usando tabu y hill-climbing, pero imponiendo las siguientes restricciones: 1. Los arcos viento --> lluvia, tormenta --> granizo y nieve --> nieveSuelo deben quedar reflejados en el DAG. 2. Ningún arco debe unir directamente la neblina con el granizo ni la niebla con la tormenta.
# Preparo la whitelist
whitelist <- matrix(c("viento", "lluvia",
"tormenta", "granizo",
"nieve", "nieveSuelo"),
ncol = 2,
dimnames = list(NULL, c("from", "to")))
# Preparo la blacklst
blacklist <- matrix(c("neblina", "granizo",
"granizo", "neblina",
"tormenta", "niebla",
"niebla", "tormenta"),
ncol = 2, byrow = TRUE,
dimnames = list(NULL, c("from", "to")))
- Vuelve a dibujar los DAG resultantes, y a partir del BIC obtenido por cada modelo, determina cuál es el mejor de todos poniéndolos en una tabla.
- Computa la fuerza de la relación entre nodos de los dos últimos modelos que combinan nuestra experiencia con el aprendizaje automático.
# Construyo el DAG "semiexperto" usando el algoritmo Tabu search
dag.tabu.semiexperto <- tabu(x = meteoro,
whitelist = whitelist,
blacklist = blacklist)
dag.tabu.semiexperto
plot(dag.tabu.semiexperto)
dag.tabu.score.semiexperto <- bnlearn::score(dag.tabu.semiexperto, data = meteoro)
dag.tabu.score.semiexperto
# Uso arc.strength para valorar la significación de los arcos
arc.strength(dag.tabu.semiexperto, data = meteoro, criterion = "x2")
# Construyo el DAG "semiexperto" usando el algoritmo hill-climbing search
dag.hc.semiexperto <- hc(x = meteoro,
whitelist = whitelist,
blacklist = blacklist)
dag.hc.semiexperto
plot(dag.hc.semiexperto)
dag.hc.score.semiexperto <- bnlearn::score(dag.hc.semiexperto, data = meteoro)
dag.hc.score.semiexperto
# Uso arc.strength para valorar la significación de los arcos
arc.strength(dag.hc.semiexperto, data = meteoro, criterion = "x2")
# Preparo una tabla de resumen con todos los resultados de los scores
summary.table <- matrix(c("DAG experto", dag.score,
"DAG Tabu", dag.tabu.score,
"DAG hill-climbing", dag.hc.score,
"DAG Tabu semi-experto", dag.tabu.score.semiexperto,
"DAG hill-climbing semi-experto", dag.hc.score.semiexperto),
ncol = 2, byrow = TRUE,
dimnames = list(NULL, c("DAG", "SCORE")))
summary.table
- Comenta brevemente los resultados obtenidos tras combinar nuestro conocimiento con el aprendizaje automático.
Como hemos visto anteriormente, el uso de algoritmos de aprendizaje ha permitido mejorar el score BIC que se había obtenido usando solo el conocimiento de un experto.
Vemos ahora que combinar las dos estrategias reduce de manera casi imperceptible la capacidad del grafo de adaptarse a los datos, pero permite evitar que relaciones probabilsticas poco significativas entren en el DAG.
Las relaciones incluidas en la whitelist:
- viento --> lluvia
- tormenta --> granizo
- nieve --> nieveSuelo
Ya aparecían en los grafos construidos a través de algoritmos de aprendizaje, que pero introducian también relaciones que hemos puesto en la blacklist, compuesta por:
- granizo <--> neblina
- tormenta <--> niebla
Como podemos ver en las tablas más abajo, la relación probabilistica entre granizo y neblina aparece, aunque poco significativa, mejorando así de manera artificial el score de los grafos, al precio de introducir mayor complejidad y cantidad de parametros. Lo mismo vale para arcos como "granizo --> rocio", que tienen una fuerza muy baja y que podemos entender como poco significativos, y que pero encontramos en los grafos construidos sin ningún tipo de conocimiento experto.
En definitiva, el método que me ha parecido más interesante ha sido el que mezcla algoritmos de aprendizaje y conocimiento previo del problema, por la siguientes razones:
- el score que obtenemos es casi tan bueno como el dado por los grafos construidos sólo con métodos de aprendizaje. La diferencia, 7 u 8 puntos sobre casi 10 mil, es básicamente nula;
- dicho score es exactamente lo mismo para los grafos con supervisión experta construidos con el algoritmo "hill-climbinig" y con el algoritmo "Tabu". Esto me parece otra prueba de la robustez del método;
- con este método, no solo hemos sido capaces de reconocer las relaciones probabilisticas más interesantes entre las caracteristicas del problema, sino que también algunas con poco sentido (e.g.: granizo --> rocio) han sido descartadas, aunque no entrasen explicitamente en la "blacklist".
arc.strength(dag.hc, data = meteoro, criterion = "x2")
arc.strength(dag.tabu, data = meteoro, criterion = "x2")